Branching process approach for Boolean bipartite networks of metabohc reactions 
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The branching process (BP) approach has been successful in explaining the avalanche dynam- 
ics in complex networks. However, its applications are mainly focused on unipartite networks, in 
which all nodes are of the same type. Here, motivated by a need to understand avalanche dynamics 
in metabolic networks, we extend the BP approach to a particular bipartite network composed of 
Boolean AND and OR logic gates. We reduce the bipartite network into a unipartite network by 
integrating out OR gates, and obtain the effective branching ratio for the remaining AND gates. 
Then the standard BP approach is applied to the reduced network, and the avalanche size distribu- 
tion is obtained. We test the BP results with simulations on the model networks and two microbial 
metabolic networks, demonstrating the usefulness of the BP approach. 

PACS numbers: 89.75.-k, 05.40.Fb, 89.20.Hh 
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Introduction — The multiplicative branching process 
(BP) approach [l[ has been successful in helping us un- 
derstand a variety of physical phenomena in complex net- 
works such as percolation [2| , epidemic spreading 0] , and 
avalanche dynamics [j|. Such an approach is valid when 
these physical phenomena do not form a nontrivial frac- 
tion of loop structures in the process of forming clusters, 
spreading diseases, and toppling cascades. Indeed, dur- 
ing the dynamic process in complex networks, the forma- 
tion of loops is a rarity, and thus, the BP approach has 
been regarded as a useful method. 

Recent studies applying the BP approach to percola- 
tion, epidemic spreading, and avalanche dynamics were 
limited to unipartite networks, in which all nodes are of 
the same type. However, there are many examples of 
bipartite networks in real-world networks, in which the 
nodes belong to one of two types. The BP approach 
has also been applied to such bipartite networks |5| and 
other multitype networks |^, 0] to study the formation 
of percolating clusters and epidemic spreading. In this 
Brief Report, we extend the previous BP formalism to 
consider a directed Boolean cascade model on a bipartite 
network composed of logic AND gates and OR gates. 
Then wc study the cascading failure problem based on 
the model, by analyzing the avalanche size distribution in 
the Boolean cascade dynamics. This generalization may 
be useful in further work for studying avalanche dynam- 
ics in bipartite or multitype networks in various systems 
and may provide a guideline for constructing a formalism 
of the BP for interacting networks [8l-[l0| . 

Boolean cascade model — The Boolean cascade model 
is defined on a directed bipartite network. The network is 
composed of two types of nodes: Boolean AND gates and 
Boolean OR gates. Each of these nodes is connected by a 
directed edge to a node of the other type, thus defining a 
bipartite network. This model was developed in previous 



studies on the basis of the reaction blockade cascade in 
metabolic networks [ll|, |l2| : When a reaction is blocked 
by, e.g., the knockout of the gcne(s) catalyzing it, its 
product cannot be produced, which in turn blocks other 
reactions, i.e., an avalanche occurs. In the dynamics, 
the metabolites act as an OR gate, since any one of the 
reactions producing that metabolite can make it turn on. 
In contrast, the reaction node is represented by the AND 
gate, since all the input metabolites should be present 
for the reaction to be activated. 

Specifically, the dynamic rule of the Boolean cascade 
model is as follows, which is also schematically illustrated 
in Fig. HJ 

(Dl) Initially, all nodes in both subsets (metabolites and 
reactions) are active, that is, their Boolean states 
are set to hi = \ for all nodes i = 1, . . . , A^. Then 
a single reaction node, say nnit, is turned off, that 

is, 6ri„it = 0- 

(D2) For all the metabolite nodes connected with 
the newly inactivated reaction nodes, update its 
Boolean state according to the Boolean OR func- 



tion, bn 



bri V 6^2 V 



Vbr 



where r,- are the 



reactions producing the metabolite m. 

(D3) For all the reaction nodes connected with the newly 
inactivated metabolite nodes, update its Boolean 
state according to the Boolean AND function, br = 



bmi/\b,n2^- 



■A6,„ 



where m, are the metabolites 
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producing the reaction r. 

(D4) Repeat (D2)-(D3) until no more inactivation oc- 
curs, and the total number of inactivated metabo- 
lite nodes, called the avalanche size, is recorded. 

(D5) Repeat the above procedure for each starting re- 
action, and obtain the avalanche size distribution 

Pais). 

Branching process approach — The propagation of the 
Boolean cascade can be understood in view of the branch- 
ing process. This approach is based on the assumption 
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FIG. 1. (Color online) Schematic illustration of avalanche 
dynamics in a Boolean bipartite network. When a reaction 
(rectangle) is deactivated, a metabolite (ellipse) would be de- 
activated if all the reactions have that metabolite as a prod- 
uct. A reaction cannot take place if any of the substrate 
metabolites is absent. Thus, metabolites and reactions cor- 
respond to logic OR and AND gates, respectively. When a 
metabolite is generated from a single reaction (with proba- 
bility pi), it can block I reactions (with probability po,m{t.)) 
and each of the metabolite branches ki {i — !,...,£) (with 
probability Po,r{ki)). The resulting value K = ki + ■ ■ ■ + ki 
corresponds to the branching number, which has probability 
qK, as derived in Eq. (1). 



that the cascade does not form a significant fraction of 
loops during propagation. This amounts to the assump- 
tion that the metabohte node should be of in-degree 1 to 
be deactivated. To start with, we assume for simplicity 
that the in-degree and out-degree of a node are uncor- 
related, that is, p{j, k) = Pi{j)po{k), where p(j, k) is the 
joint probability distribution of the in-degree j and the 
out-degree k. 

Next, we reduce the bipartite network to a unipartite 
network composed only of reactions. In this transfor- 
mation, the probability distribution qK of the effective 
branching ratio K for a given reaction is derived. This 
can be more easily obtained in terms of the generating 
function, 



Q(w) = Y, IK^'' 



A'=0 



oo oo 



(1-Pl)+Pl^ ^ Po,m(^) JJpo,r(fc,)w^'=i 
e=0 ki,--- M=0 i=l 



= (1 -pi) + Pl^Po.rnWQriujf 



e=o 



(1 -pi) +PlQrn{Qr{^)), 



(1) 



where pi is the probability that the in-degree of a 
metabolite encountered is one (see Fig. 1), and the fc^'s 
and i are randomly chosen from the probability distribu- 
tions po,r{k) and Po,m,ik), respectively. Then, the effec- 
tive branching ratio K is determined by K ~ ki -\- k2 + 
■ ■ ■ + ki. Qm{'^) and Qr{i-^) are the generating functions 



defined as 

oo oo 

Q„i{w) ^ Ypo,mik)u;'' and (9r(w) = ^Po,r(fc) 



U! 



fc=0 



fe=0 



(2) 
The criticality condition is (K) = Q'(w)|a;=i = 1. 

Next, we solve for the avalanche size distribution Pais) 
of the Boolean dynamics with the assumption that it pro- 
ceeds in a branching tree pattern. Using the standard 
steps of the BP approach, we write the avalanche size 
distribution Pa (s) as 



Pais) = 



qK 



K=0 si = l 



51 ■•■ 51 Paisi) ■ ■ ■ Pa{sK)Sj2^^^s,^s-V 



(3) 

where a seed deactivated reaction has K branches, and 
from each of the branches, s^-sized successive deactivated 
reactions follow under the condition ^^ Si ~ s — I. This 
formula can be rewritten in terms of the generating func- 
tions Pa{z) = J2T=iPa{s)z'' and Q{uj) as 



Va{z) = zQiVaiz)). 



(4) 



Assuming that the avalanche size distribution follows a 
power law, i.e., Pai-s) ~ s~'^ , we can expand the generat- 
ing function near 2 = 1 as 



Paiz) C^l-{l-z)+Ci{l-zY-'+C2{l-zf + - 



(5) 



Let us consider the criticality condition Q'{id)\uj=i = 1. 
Applying it to Eq. (1), we obtain the mean branching 
ratio as 



Q'i^)l 



PlQ'„^{Qr{l))Q'Ai) ^ Pl{ko,m){ko.r}, (6) 



where {ko.m) and {ko.r) arc the mean out-degrees of the 
metabolites and the reactions, respectively. Therefore, 
we have to tune the parameter pi to maintain the criti- 
cality condition, pi is the probability of reaching a node 
(metabolite) whose inward degree is one, which is given 
by the formula pk = kpi^„iik)/{ki^rn) with fc = 1 as 
Pi = Pi^m{^) / {ki.m) , so that the mean branching ratio 
becomes 



{K)^ 






(7) 



Power-law degree distributions — As an example, we 
consider the cases where Po,m{k) and po,r{k) follow power 
laws, i.e., Po,m{k) ~ k~''"^ and Po,r(^) ^ k~'^'' . As shown 
below, this is an interesting regime not only because we 
have nontrivial scaling behavior in Pais), but also being 
relevant for real- world metabolic networks [l3| . 

In this case, we can determine the expansion of Q{uj) 
near a; sa 1 from the smaller of 7„i and 7^, i.e., 7 = 
min[7„i,7r], as 



Q(w) ~ l-(l-cj) 



Al(l-w)''-l 


(2 < 7 < 3) 


-A2il-ujfHl- 


-uj) (7 = 3), 


^3(1-^)2 


(7 > 3). 




(8) 



To solve the self-consistent equation ([4]), we set lo 
Va{z), and consequently, z = V^^ = u/Q{lo) « 1 
Ai{l - w)^-'^ H . Thus, we obtain 



7/(7- 
3/2 



1) (2<7<3), 
(7 > 3). 



(9) 



For the case in which all pi^„i(k), po^mik), and Po,r{k) 
follow the same power-law distribution with the same 
exponent 7, (A') = C(7— 1)/C(7)^> which is larger than 1 
for finite 7. Thus, the BP shows a supercritical behavior. 
In the supercritical regime (A') > 1 Jl4l], the avalanche 
size distribution behaves as 



Pa{s) 






(10) 



in which Sc ^ \ (AT) — 1 1 
2 < 7 < 3, and 

Pa{s) 



with a = (7 - l)/(7 - 2) for 



5-3/2g-./. 



(11) 



withsc-^ |(A:)-1|-2 for 7 > 3. 

Numerical simulations — To test the BP predictions 
with simulation results, we first construct a model 
Boolean bipartite network, which is a generalization of 
configuration model, as follows: 

(51) The system is composed of A^ nodes, which are 
composed of two subsets, Ad and 7?,, of equal size, 

N/2. 

(52) For a node v € AA, assign in-degree ki^y and out- 
degree ko^v randomly from the probability distri- 
butions pi^m{k) and po.m{k), respectively. Simi- 
larly, for a node w ^ TZ, assign in-degree fci_„, and 
out-degree ko,w from the probability distributions 
Pi^r{k) and Po,rik), respectively. 

(53) Choose an ordered pair of vertices {v, w) and join 
them by the directed edge w — )• w if they belong to 
different subsets, provided that there remain a free 
outward arrow at v and a free inward arrow at w 
and that the pair is not connected yet, that is, we 
disallow multiple arrows between the node pair. 

(54) Repeat (S3) until there is no stub left. 

(55) Assign Boolean OR gates to all nodes in the subset 
Ai and assign Boolean AND gates to those in TZ. 

Simulation results for the Boolean cascade dynamics 
on the model networks with in- and out-degree distribu- 
tion of same power-law form, with 7 = 2.5, 3.5, and 4.5 
are shown in Fig. 2, showing a good agreement with the 
BP prediction. 

Real-world metabolic networks — We now consider the 
Boolean cascade dynamics on the real-world metabolic 
networks of Escherichia coli and Saccharomyces cere- 
visiae [1^, which contain 1,188 metabolites and 1,489 




FIG. 2. (Color online) The avalanche size distribution Pa{s) 
for tlie Boolean cascade model on the directed bipartite 
scale-free network with the in- and out-degree distribution 
Pdik) ~ k ' and 7 = 2.5 (red), 3.5 (green), and 4.5 (bhie). 
The curves are a good fit to the theoretical function Pa{s) = 
yls~^exp(— s/sc) with r values of 1.71 (red), 1.56 (green), 
and 1.47 (blue). Note that we are in the supercritical regime 
so that isolated peaks at large s appear. The system size is 
iV = 10^ 



reactions for E. coli and 680 metabolites and 852 reac- 
tions for S. cerevisiae, respectively. The out-degree dis- 
tributions for the metabolites and the reactions in the 
metabolic networks are measured. For both networks, 
the out-degree distribution for metabolites is found to 
follow approximately a power law with the exponent of 
7o,m ~ 2.2 (Fig. 3, insets), in agreement with previous 
results [13. On the other hand, the out-degree distribu- 
tion for reactions is found to decay exponentially (Fig. 3, 
insets). 

We run the Boolean cascade dynamics on these net- 
works and measure the avalanche size distribution Pai-s)- 
From numerical simulation, we found that the avalanche 
size distribution for both species has an approximate 
power-law tail, with the exponent r « 1.8 (Fig. 3). This 
power-law exponent is not far from the BP prediction by 
Eq. (9) with 7 = 2.2. 

To further compare the numerical simulation results 
with BP predictions, we numerically solved Eq. (4) by 
plugging in the empirical degree distributions, the result 
of which is depicted also in Fig. 3. The theoretical Pa{s) 
deviates from a pure power law, because the dynamics is 
super-critical, as (AT) w 3.26 (A. coli) and 3.85 (5. cere- 
visiae), respectively. Also, real- world metabolic networks 
contain the degree-degree correlation and clustering, as 
well as are finite, which lead to deviations from the mean- 
field type analytic predictions [16j . For example, we note 
the peaks in the empirical Pa(s) near s « 20 ~ 30, that 
are absent in theoretical Pa{s), which might arise due 
to the existence of compact modules or cycles in the 
metabolic networks, such as TCA cycle. Despite these 
complications, the power-law behaviors obtained from 
numerical simulations and BP calculations agree reason- 




ably, demonstrating the usefulness of BP approach to 
understanding the scaling behaviors. 



Summary — In conclusion, we have studied the 
Boolean cascade dynamics occurring in directed bi- 
partite networks of Boolean AND and OR logic gates, 
inspired by the cascading failure of reactions in metabolic 
networks. A branching process approach was developed 
to study this type of dynamics. Theoretically obtained 
exponent for the avalanche size distribution is in good 
agreement with the simulation results for the model 
network as well as for two empirical microbial metabolic 
networks. 



FIG. 3. (Color online) The avalanche size distributions on 
real-world metabolic networks, (a) is for E. coli and (b) is for 
S. cerevisiae. Circles are for simulation results, and squares 
are for BP predictions with empirical degree distributions. 
Solid lines are guidelines with the slope of —1.83 which sat- 
isfies the relation r = 7/(7 — 1) with 7 = 2.2. Insets are 
out-degree distributions of metabolites (triangles) and reac- 
tions (diamonds). Solid guidelines in the insets have a slope 
-2.2. 
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